---
title: "Simonsen_2020_reproduction"
author: "S. Shin"
date: "2/2/2022"
output: html_document
---

---
title: "Distinguishing Oneself from Others - Analysis Script"
author: "Arndis Simonsen & Riccardo Fusaroli"
date: "9/8/2019"
output: html_document
---

# Analysis script for the Distinguisting Oneself from Others manuscript

## Structure of the markdown

* Load libraries; load data and preprocess it
* Descriptive statistics
* Modeling of group differences and hypothesis testing
* Modeling of individual differences: TASIT
* Modeling of individual differences: TASIT and IQ
* Modeling of individual differences: Animated Triangles (ATT)
* Modeling of individual differences: Clinical features
* Modeling of individual differences: Clinical features and TASIT
* Figure 4: Counterfactual plot of relation between altercentric intruision, TASIT, ATT and relevant symptoms


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)
set.seed(2020)

## Function to  reduce the amount of code to fit models
brmsFun<-function(Name, formula, data,family = gaussian, prior = SkepticPrior1, iter = 4000, path="models"){
  m <- brm(formula, 
           data = data,
           family = family,
           iter = iter, 
           chains = 2, 
           cores = 2,
           prior = prior,
           sample_prior = T,
           file = here(path, Name),
           control = list(adapt_delta = 0.99,
                         max_treedepth = 20))

  m <- add_criterion(m,criterion=c("waic","loo", "bayes_R2"),reloo=F,file = here(path, Name))
  return(m)
}

```

## Load data and libraries


```{r load data and libraries}

pacman::p_load(tidyverse,
               see,
               brms,
               here,
               bayestestR,
               patchwork,
               viridis,
               ggridges,
               patchwork)

dT <- read_csv("Simonsen_2020.csv") %>% mutate(
  Consistency = as.factor(Consistency),
  Group = as.factor(Group),
  Self = as.factor(Self),
  ATT = as.ordered(ATT)
)


## Remove practice trials
dT <- dT %>%
  subset(Trial > 9 & True == "j") # Remove practice trials and correct no responses

## Calculate how many rt's below 200 ms in the 2 groups
dT %>% group_by(Group) %>% summarize(n())
dT %>% group_by(Group) %>% summarize(sum(RT>0 & RT<200))
dT %>% group_by(Group) %>% summarize(sum(RT==0))

cat(paste0(
  "Probability of timeouts in controls = ", round(3/828, 4), "\n",
  "Probability of timeouts in schizophrenia = ", round(32/792, 4), "\n",
  "Probability of RT below 200 (in schizophrenia) = ", round(1/1620, 4)
))

## Excluding all RTs <200 ms. 
dT <- subset(dT,RT > 200)

## Calculate errors
dT %>% group_by(Group) %>% summarize(n())
dT %>% group_by(Group) %>% summarize(sum(Accuracy == 0))

cat(paste0(
  "Probability of errors in controls = ", round(28/828, 4), "\n",
  "Probability of errors in schizophrenia = ", round(80/792, 4)
))


## Double check RTs
xx <- dT %>% group_by(partID,Subject,Group,MatchID,Self,Consistency) %>% dplyr::summarise(RT = mean(RT,na.rm=T))
xx

## Check how many trials by participant
dT %>% group_by(Group, partID) %>% summarize(TrialCount=n()) %>% group_by(Group) %>% summarize(mean(TrialCount))

```

## Descriptive statistics

```{r Descriptive statistics}
## Summing up participants' features by group

x<-dT %>% 
  group_by(Subject, Group) %>% 
  dplyr::summarise(Relevant = mean(Relevant, na.rm = T), 
                   Irrelevant = mean(Irrelevant, na.rm = T),
                   Tasit = mean(Tasit, na.rm=T),
                   ATT = mean(as.numeric(ATT), na.rm=T),
                   EduYear = mean(EduYear, na.rm=T),
                   IQ = mean(EsToIQ, na.rm=T),
                   DuraOPUS = mean(DuraOPUS)) %>%
  group_by(Group) %>%
  dplyr::summarise(RelevantM = mean(Relevant, na.rm=T),
                   Relevant25 = quantile(Relevant, 0.025, na.rm=T),
                   Relevant975 = quantile(Relevant, 0.975, na.rm=T),
                   IrrelevantM = mean(Irrelevant, na.rm=T),
                   Irrelevant25 = quantile(Irrelevant, 0.025, na.rm=T),
                   Irrelevant975 = quantile(Irrelevant, 0.975, na.rm=T),
                   TasitM = mean(Tasit, na.rm=T),
                   Tasit25 = quantile(Tasit, 0.025, na.rm=T),
                   Tasit975 = quantile(Tasit, 0.975, na.rm=T),
                   TrianglesM = mean(as.numeric(ATT), na.rm=T),
                   Triangles25 = quantile(as.numeric(ATT), 0.025, na.rm=T),
                   Triangles975 = quantile(as.numeric(ATT), 0.975, na.rm=T),
                   EduYearM= mean(EduYear, na.rm=T),
                   Duration = mean(DuraOPUS))
x
## Test for differences in Education and IQ
TestD <- dT %>% 
  group_by(Subject, Group) %>% 
  dplyr::summarise(Relevant = mean(Relevant, na.rm = T), 
                   Irrelevant = mean(Irrelevant, na.rm = T),
                   Tasit = mean(Tasit, na.rm=T),
                   ATT = mean(as.numeric(ATT), na.rm=T),
                   EduYear = mean(EduYear, na.rm=T),
                   IQ = mean(EsToIQ, na.rm=T),
                   Age = mean(Age, na.rm=T),
                   Gender = Gender[1])

EduIQ <- brm(EduYear ~ mo(IQ),
             family = cumulative,
             subset(TestD, Group=="Schizophrenia"),
             chain = 2)


Edu_m <- brm(EduYear ~ 0 + Group,
             family = poisson,
             TestD,
             prior = prior(normal(2.6, 0.2), class = b),
             sample_prior = T,
             chains = 2)
pp_check(Edu_m, nsamples=100)

hypothesis(Edu_m, "GroupControl = GroupSchizophrenia")
hypothesis(Edu_m, "GroupControl > GroupSchizophrenia")

IQ_m <- brm(IQ ~ 0 + Group,
             family = poisson,
             TestD,
             prior = prior(normal(4.6, 0.25), class = b),
             sample_prior = T,
             chains = 2)
pp_check(IQ_m, nsamples=100)

hypothesis(IQ_m, "GroupControl = GroupSchizophrenia")
hypothesis(IQ_m, "GroupControl > GroupSchizophrenia")

Tasit_m <- brm(Tasit ~ 0 + Group,
             family = poisson,
             TestD,
             prior = prior(normal(4.2, 0.18), class = b),
             sample_prior = T,
             chains = 2)
pp_check(Tasit_m, nsamples=100)

hypothesis(Tasit_m, "GroupControl = GroupSchizophrenia")
hypothesis(Tasit_m, "GroupControl > GroupSchizophrenia")

ATT_m <- brm(ATT ~ 0 + Group,
             family = poisson,
             TestD,
             prior = prior(normal(2.43, 0.6), class = b),
             sample_prior = T,
             chains = 2)
pp_check(ATT_m, nsamples=100)

hypothesis(ATT_m, "GroupControl = GroupSchizophrenia")
hypothesis(ATT_m, "GroupControl > GroupSchizophrenia")
```


## Statistical modeling of the group differences

```{r Statistical modeling of the group differences}

## Define model formulas
GroupAccuracy_f <- bf(
  Accuracy ~ 0 + Self : Consistency : Group + (0 + Self : Consistency : Group | MatchID) + (0 +  Group |Picture)
)

GroupRT_f <- bf(
  RT ~ 0 + Self : Consistency : Group + (0 + Self : Consistency : Group | MatchID) + (0 +  Group |Picture)
)

## Define priors
get_prior(GroupAccuracy_f,
          data = dT,
          family = "binomial")

priorAccuracy <- c(
  prior(normal(1,1), class=b),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

get_prior(GroupRT_f,
          data = dT,
          family = shifted_lognormal)

priorRT <- c(
  prior(normal(6,.5), class=b),
  prior(lkj(3), class=cor),
  prior(normal(0,.1), class=sd),
  prior(normal(0,.3), class=sigma)
)

## Prior predictive checks

GroupAccuracy_m1_PriorCheck <- brm(
  GroupAccuracy_f,
  data = dT,
  family = bernoulli,
  prior = priorAccuracy,
  sample_prior = "only",
  file = here("models","GroupAccuracy_m1_PriorCheck"),
  cores = 2,
  chains = 2,
  iter = 4000,
  control = list(
    max_treedepth = 20,
    adapt_delta = 0.99)
)

GroupRT_m1_PriorCheck <- brm(
  GroupRT_f,
  data = subset(dT, Accuracy == "1"),
  family = shifted_lognormal,
  prior = priorRT,
  sample_prior = "only",
  file = here("models","GroupRT_m1_PriorCheck"),
  cores = 2,
  chains = 2,
  iter = 4000,
  control = list(
    max_treedepth = 20,
    adapt_delta=0.99)
)

## Plot prior predictive checks
pp_check(GroupAccuracy_m1_PriorCheck, nsamples = 100)
pp_check(GroupRT_m1_PriorCheck, nsamples = 100)

## Fit model
GroupAccuracy_m1 <- brmsFun(
  Name = "GroupAccuracy_m1", 
  formula = GroupAccuracy_f, 
  data = dT,
  family = bernoulli,
  prior = priorAccuracy)

## Posterior predictive checks
pp_check(GroupAccuracy_m1, nsamples = 100)

## Check model diagnostics
summary(GroupAccuracy_m1)

GroupRT_m1 <- brmsFun(
  Name = "GroupRT_m1", 
  formula = GroupRT_f, 
  data = subset(dT,Accuracy == "1"),
  family = shifted_lognormal,
  prior = priorRT)

## Posterior predictive checks
pp_check(GroupRT_m1, nsamples = 100)

## Check model diagnostics
summary(GroupRT_m1)

## Sanity check of main effects

hypothesis(GroupAccuracy_m1, c(
  "(SelfOther:Consistency0:GroupControl + SelfOther:Consistency1:GroupControl + SelfOther:Consistency0:GroupSchizophrenia + SelfOther:Consistency1:GroupSchizophrenia)/4 > (SelfSelf:Consistency0:GroupControl + SelfSelf:Consistency1:GroupControl + SelfSelf:Consistency0:GroupSchizophrenia + SelfSelf:Consistency1:GroupSchizophrenia)/4",
  "(SelfOther:Consistency0:GroupControl + SelfOther:Consistency1:GroupControl + SelfOther:Consistency0:GroupSchizophrenia + SelfOther:Consistency1:GroupSchizophrenia)/4 > 0",
  "(SelfSelf:Consistency0:GroupControl + SelfSelf:Consistency1:GroupControl + SelfSelf:Consistency0:GroupSchizophrenia + SelfSelf:Consistency1:GroupSchizophrenia)/4>0"
))

hypothesis(GroupRT_m1, c(
  "(SelfOther:Consistency0:GroupControl + SelfOther:Consistency1:GroupControl + SelfOther:Consistency0:GroupSchizophrenia + SelfOther:Consistency1:GroupSchizophrenia)/4 < (SelfSelf:Consistency0:GroupControl + SelfSelf:Consistency1:GroupControl + SelfSelf:Consistency0:GroupSchizophrenia + SelfSelf:Consistency1:GroupSchizophrenia)/4",
  "(SelfOther:Consistency0:GroupControl + SelfOther:Consistency1:GroupControl + SelfOther:Consistency0:GroupSchizophrenia + SelfOther:Consistency1:GroupSchizophrenia)/4 > 0",
  "(SelfSelf:Consistency0:GroupControl + SelfSelf:Consistency1:GroupControl + SelfSelf:Consistency0:GroupSchizophrenia + SelfSelf:Consistency1:GroupSchizophrenia)/4>0"
))

```

## Hypothesis testing
```{r Hypothesis testing for group differences}

print("Testing for egocentric intrusion")

x1 <- hypothesis(GroupAccuracy_m1, "SelfOther:Consistency0:GroupControl < SelfOther:Consistency1:GroupControl")
x2 <- hypothesis(GroupRT_m1, "SelfOther:Consistency0:GroupControl > SelfOther:Consistency1:GroupControl")
x3 <- hypothesis(GroupAccuracy_m1, "SelfOther:Consistency1:GroupSchizophrenia - SelfOther:Consistency0:GroupSchizophrenia  <
                 SelfOther:Consistency1:GroupControl - SelfOther:Consistency0:GroupControl")
x4 <- hypothesis(GroupRT_m1, "SelfOther:Consistency1:GroupSchizophrenia - SelfOther:Consistency0:GroupSchizophrenia  <
                 SelfOther:Consistency1:GroupControl - SelfOther:Consistency0:GroupControl")

text <- paste0("As expected, when healthy controls had to take the avatar’s perspective, their own perspective interfered, i.e. they were slower and made more errors on inconsistent trials (egocentric intrusion - Accuracy: ",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), "; RT: ",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ". ",
              "However, this was also the case for patients and the two groups were not credibly different (Group difference in egocentric intrusion for Accuracy: ",
               "β = ", round(x3$hypothesis$Estimate, 2),
              ", SE = ", round(x3$hypothesis$Est.Error, 2),
              ", CI (", round(x3$hypothesis$CI.Lower, 2), ", ", round(x3$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x3$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x3$hypothesis$Post.Prob, 2), "; RT: ",
              "β = ", round(x4$hypothesis$Estimate, 2),
              ", SE = ", round(x4$hypothesis$Est.Error, 2),
              ", CI (", round(x4$hypothesis$CI.Lower, 2), ", ", round(x4$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x4$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x4$hypothesis$Post.Prob, 2), ").")

cat(text)

print("Testing for altercentric intrusion")

x1 <- hypothesis(GroupAccuracy_m1, "SelfSelf:Consistency0:GroupControl < SelfSelf:Consistency1:GroupControl")
x2 <- hypothesis(GroupRT_m1, "SelfSelf:Consistency0:GroupControl < SelfSelf:Consistency1:GroupControl")
x3 <- hypothesis(GroupAccuracy_m1, "SelfSelf:Consistency1:GroupSchizophrenia - SelfSelf:Consistency0:GroupSchizophrenia  >
                 SelfSelf:Consistency1:GroupControl - SelfSelf:Consistency0:GroupControl")
x4 <- hypothesis(GroupRT_m1, "SelfSelf:Consistency1:GroupSchizophrenia - SelfSelf:Consistency0:GroupSchizophrenia  <
                 SelfSelf:Consistency1:GroupControl - SelfSelf:Consistency0:GroupControl")

text <- paste0("When the controls had to take their own perspective, the avatar’s perspective also interfered with task performance (altercentric intrusion). Specifically, controls made more errors when the avatar’s perspective was inconsistent with their own. However, the RTs were not different in the two conditions (Accuracy: ",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), "; RT: ",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ". ",
              "Patients on the other hand, displayed more altercentric intrusion compared to controls, as reflected in the higher number of errors and slower RTs when the two perspectives differed (Accuracy: ",
               "β = ", round(x3$hypothesis$Estimate, 2),
              ", SE = ", round(x3$hypothesis$Est.Error, 2),
              ", CI (", round(x3$hypothesis$CI.Lower, 2), ", ", round(x3$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x3$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x3$hypothesis$Post.Prob, 2), "; RT: ",
              "β = ", round(x4$hypothesis$Estimate, 2),
              ", SE = ", round(x4$hypothesis$Est.Error, 2),
              ", CI (", round(x4$hypothesis$CI.Lower, 2), ", ", round(x4$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x4$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x4$hypothesis$Post.Prob, 2), ").")

cat(text)

# Generate posteriors for plots

PosteriorAccuracy <- posterior_samples(
  GroupAccuracy_m1, 
  pars = c(
  "b_SelfOther:Consistency0:GroupControl",
  "b_SelfSelf:Consistency0:GroupControl",
  "b_SelfOther:Consistency1:GroupControl",
  "b_SelfSelf:Consistency1:GroupControl",
  "b_SelfOther:Consistency0:GroupSchizophrenia",
  "b_SelfSelf:Consistency0:GroupSchizophrenia",
  "b_SelfOther:Consistency1:GroupSchizophrenia",
  "b_SelfSelf:Consistency1:GroupSchizophrenia")
  ) %>%
  gather(condition, posterior, 1:8) %>%
  mutate(
    Self = ifelse(grepl("SelfSelf",condition),"Self","Other"),
    Consistency = factor(ifelse(grepl("Consistency0",condition),"Inconsistent","Consistent"), levels = c("Inconsistent", "Consistent")),
    Group = ifelse(grepl("GroupControl",condition),"Control","Schizophrenia"),
    Accuracy = inv_logit_scaled(posterior),
    condition = NULL
  )

PosteriorRT <- posterior_samples(
  GroupRT_m1, 
  pars = c(
  "b_SelfOther:Consistency0:GroupControl",
  "b_SelfSelf:Consistency0:GroupControl",
  "b_SelfOther:Consistency1:GroupControl",
  "b_SelfSelf:Consistency1:GroupControl",
  "b_SelfOther:Consistency0:GroupSchizophrenia",
  "b_SelfSelf:Consistency0:GroupSchizophrenia",
  "b_SelfOther:Consistency1:GroupSchizophrenia",
  "b_SelfSelf:Consistency1:GroupSchizophrenia",
  "ndt")
  ) %>%
  gather(condition, posterior, 1:8) %>%
  mutate(
    lndt = log(ndt),
    RT = exp(posterior) + ndt,
    Self = ifelse(grepl("SelfSelf",condition),"Self","Other"),
    Consistency = factor(ifelse(grepl("Consistency0",condition),"Inconsistent","Consistent"), levels = c("Inconsistent", "Consistent")),
    Group = ifelse(grepl("GroupControl",condition),"Control","Schizophrenia"),
    condition = NULL,
    prior_ndt = NULL
  )

x1 <- data.frame(
  Group = rep(c("Control","Schizophrenia"), each = 4000),
  AltercentricIntrusionAccuracy = c(
    subset(PosteriorAccuracy, Group=="Control" & Self=="Self" & Consistency=="Consistent")$Accuracy - 
     subset(PosteriorAccuracy, Group=="Control" & Self=="Self" & Consistency=="Inconsistent")$Accuracy,
    subset(PosteriorAccuracy, Group=="Schizophrenia" & Self=="Self" & Consistency=="Consistent")$Accuracy - 
     subset(PosteriorAccuracy, Group=="Schizophrenia" & Self=="Self" & Consistency=="Inconsistent")$Accuracy),
  EgocentricIntrusionAccuracy = c(
    subset(PosteriorAccuracy, Group=="Control" & Self=="Other" & Consistency=="Consistent")$Accuracy - 
     subset(PosteriorAccuracy, Group=="Control" & Self=="Other" & Consistency=="Inconsistent")$Accuracy,
    subset(PosteriorAccuracy, Group=="Schizophrenia" & Self=="Other" & Consistency=="Consistent")$Accuracy - 
     subset(PosteriorAccuracy, Group=="Schizophrenia" & Self=="Other" & Consistency=="Inconsistent")$Accuracy),
  AltercentricIntrusionRT = c(
    subset(PosteriorRT, Group=="Control" & Self=="Self" & Consistency=="Inconsistent")$RT - 
     subset(PosteriorRT, Group=="Control" & Self=="Self" & Consistency=="Consistent")$RT,
    subset(PosteriorRT, Group=="Schizophrenia" & Self=="Self" & Consistency=="Inconsistent")$RT - 
     subset(PosteriorRT, Group=="Schizophrenia" & Self=="Self" & Consistency=="Consistent")$RT),
  EgocentricIntrusionRT = c(
    subset(PosteriorRT, Group=="Control" & Self=="Other" & Consistency=="Inconsistent")$RT - 
     subset(PosteriorRT, Group=="Control" & Self=="Other" & Consistency=="Consistent")$RT,
    subset(PosteriorRT, Group=="Schizophrenia" & Self=="Other" & Consistency=="Inconsistent")$RT - 
     subset(PosteriorRT, Group=="Schizophrenia" & Self=="Other" & Consistency=="Consistent")$RT)
)


AlloAccu_p1 <- ggplot(
  subset(PosteriorAccuracy, Self=="Self"), aes(x = Consistency, y = Accuracy,  color=Group)) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult=1), 
                 geom = "pointrange", position = position_dodge(width = 1), show.legend = FALSE) + 
  theme_classic() +
  xlab("") + 
  ggtitle('Accuracy')  +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_brewer(palette="Dark2")

AlloRT_p1 <- ggplot(
  subset(PosteriorRT,Self=="Self"), aes(x = Consistency, y = RT,  color=Group)) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult=1), 
                 geom = "pointrange", position = position_dodge(width = 1)) + 
  theme_classic() +
  xlab("") +
  ylab("Reaction time (ms)") +
  ggtitle('RT')  +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_brewer(palette="Dark2")


AlloAccu_p2 <- ggplot(x1, aes(x = AltercentricIntrusionAccuracy,  color=Group, fill=Group, group=Group)) +
  geom_density(alpha=0.3, show.legend = FALSE) +
  theme_classic() +
  ylab("") +
  xlab("") +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme(axis.text.y = element_blank())

AlloRT_p2 <- ggplot(x1, aes(x = AltercentricIntrusionRT,  color=Group, fill=Group, group=Group)) +
  geom_density(alpha=0.3) +
  theme_classic() +
  xlab("") +
  ylab("") +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme(axis.text.y = element_blank())

EgoAccu_p1 <- ggplot(
  subset(PosteriorAccuracy, Self=="Other"), aes(x = Consistency, y = Accuracy,  color=Group)) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult=1), 
                 geom = "pointrange", position = position_dodge(width = 1), show.legend = FALSE) + 
  theme_classic() +
  xlab("") +
  ggtitle('Accuracy')  +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_brewer(palette="Dark2")

EgoRT_p1 <- ggplot(
  subset(PosteriorRT,Self=="Other"), aes(x = Consistency, y = RT,  color=Group)) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult=1), 
                 geom = "pointrange", position = position_dodge(width = 1), show.legend = FALSE) + 
  theme_classic() +
  xlab("") +
  ylab("Reaction time (ms)")+
  ggtitle('RT') +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_color_brewer(palette="Dark2")

EgoAccu_p2 <- ggplot(x1, aes(x = EgocentricIntrusionAccuracy,  color=Group, fill=Group, group=Group)) +
  geom_density(alpha=0.3, show.legend = FALSE) +
  theme_classic() +
  ylab("") +
  xlab("") +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme(axis.text.y = element_blank())

EgoRT_p2 <- ggplot(x1, aes(x = EgocentricIntrusionRT,  color=Group, fill=Group, group=Group)) +
  geom_density(alpha=0.3, show.legend = FALSE) +
  theme_classic() +
  xlab("") +
  ylab("") +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme(axis.text.y = element_blank())

EgoAccu_p1 + EgoRT_p1 +  AlloAccu_p1 + AlloRT_p1 +
  EgoAccu_p2 + EgoRT_p2 + AlloAccu_p2 + AlloRT_p2 +
  plot_layout(nrow = 2, heights=c(3,1)) + 
  plot_annotation(title = 'Egocentric Intrusion                                                                                 Altercentric Intrusion') +
  plot_annotation(tag_levels = 'A') & theme(plot.title = element_text(hjust = 0.5),  
  plot.tag = element_text(size = 8))

ggsave(here("plots","EgocentricIntrusion.svg"), width = 16, height = 6)

# Spaghetti plots??
       
dT <- dT %>% mutate(
  ConsistencyN = as.numeric(Consistency)-1,
  AccuracyN = as.numeric(Accuracy),
  MatchID = as.factor(MatchID))

ggplot(dT,aes(ConsistencyN, AccuracyN, group = MatchID, color = MatchID)) +
  geom_smooth(method=lm, se=F, alpha=0.3) +
  facet_wrap(Self ~ Group) +
  theme_classic()

ggsave(here("plots","Accuracy_spaghetti.svg"))

ggplot(subset(dT,AccuracyN==1),aes(ConsistencyN, RT, group = MatchID, color = MatchID)) +
  geom_smooth(method=lm, se=F) +
  facet_wrap(Self ~ Group) +
  theme_classic()
ggsave(here("plots","RT_spaghetti.svg"))

hypothesis(GroupAccuracy_m1, "SelfSelf:Consistency1:GroupSchizophrenia - SelfSelf:Consistency0:GroupSchizophrenia  >
                 SelfSelf:Consistency1:GroupControl - SelfSelf:Consistency0:GroupControl", group = "MatchID", scope = "coef")


```

## Investigating individual differences in altercentric intrusion 

* TASIT
* ATT
* Relevant psychotic symptoms

### TASIT

```{r Tasit}
TasitAccuracy_mo_f <- bf(
  Accuracy ~ 0 + Self : Consistency + mo(TotAB) : Self : Consistency + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)
TasitRT_mo_f <- bf(
  RT ~ 0 + Self : Consistency + mo(TotAB) : Self : Consistency + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorTasitAccuracy_mo <- c(
  prior(normal(1, 1), class=b),
  prior(normal(0, .3), class=b, coef=SelfOther:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moTotAB),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

priorTasitRT_mo <- c(
  prior(normal(6,.5), class=b),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moTotAB),
  prior(lkj(3), class=cor),
  prior(normal(0,.1), class=sd),
  prior(normal(0,.3), class=sigma)
)

TasitAccuracy_mo_m1 <- brmsFun(
  Name = "TasitAccuracy_mo_m1",
  formula = TasitAccuracy_mo_f,
  data = subset(dT, Group=="Schizophrenia"),
  family = bernoulli,
  prior = priorTasitAccuracy_mo
)
TasitAccuracy_mo_m1_ctrl <- brmsFun(
  Name = "TasitAccuracy_mo_m1_ctrl",
  formula = TasitAccuracy_mo_f,
  data = subset(dT, Group=="Control"),
  family = bernoulli,
  prior = priorTasitAccuracy_mo
)

# Is there an interaction?
hypothesis(TasitAccuracy_mo_m1, "80 *  SelfSelf:Consistency1:moTotAB - 80 * SelfSelf:Consistency0:moTotAB  > 0", class="bsp")
hypothesis(TasitAccuracy_mo_m1_ctrl, "80 *  SelfSelf:Consistency1:moTotAB - 80 * SelfSelf:Consistency0:moTotAB  > 0", class="bsp")

# Calculating allocentric effect with random intentional score of 0
hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency0 > 0")
hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency1 > 0")
inv_logit_scaled(1.36) - inv_logit_scaled(1.32)

hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency0:moTotAB > 0", class="bsp")
hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency1:moTotAB > 0", class="bsp")
inv_logit_scaled(1.36 + 40 * 0.07) - inv_logit_scaled(1.32 + 40 * 0.01)
inv_logit_scaled(1.36 + 80 * 0.07) - inv_logit_scaled(1.32 + 80 * 0.01)

TasitRT_mo_m1 <- brmsFun(
  Name = "TasitRT_mo_m1",
  formula = TasitRT_mo_f,
  data = subset(dT,Accuracy == "1" & Group=="Schizophrenia"),
  family = shifted_lognormal,
  prior = priorTasitRT_mo
)

Preds <- predict(TasitAccuracy_mo_m1)
dT$Preds <- NULL
dT$Preds[dT$Group=="Schizophrenia"] <- Preds[,1]

ggplot(dT, aes(TotAB,Preds, color=Consistency)) + 
  geom_point() + 
  geom_smooth(method=lm, formula=y~poly(x,2)) +
  facet_wrap(.~Self) +
  theme_classic()

PredsRT <- predict(TasitRT_mo_m1)
dT$PredsRT <- NULL
dT$PredsRT[dT$Group=="Schizophrenia" & dT$Accuracy==1] <- PredsRT[,1]

ggplot(dT, aes(TotAB,PredsRT, color=Consistency)) + 
  geom_point() + 
  geom_smooth(method=lm, formula=y~poly(x,2)) +
  facet_wrap(.~Self) +
  theme_classic()

# Is there an interaction?
hypothesis(TasitRT_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  < 0", class="bsp")


###
# 

x1 <- hypothesis(TasitAccuracy_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  < 0", class="bsp")
x2 <- hypothesis(TasitRT_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  < 0", class="bsp")

x3a <- hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency0 > 0")
x3aa <- hypothesis(TasitAccuracy_mo_m1, "40 * SelfSelf:Consistency0:moTotAB > 0", class="bsp")
x3b <- hypothesis(TasitAccuracy_mo_m1, "80 * SelfSelf:Consistency0:moTotAB > 0", class="bsp")
x4a <- hypothesis(TasitAccuracy_mo_m1, "SelfSelf:Consistency1 > 0")
x4aa <- hypothesis(TasitAccuracy_mo_m1, "40 * SelfSelf:Consistency1:moTotAB > 0", class="bsp")
x4b <- hypothesis(TasitAccuracy_mo_m1, "80 * SelfSelf:Consistency1:moTotAB > 0", class="bsp")

text <- paste0("TASIT total score was credibly related to the altercentric intrusion effect for accuracy (",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), "), but not RT performance (",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ".",
              "In particular, the higher the TASIT score, the larger the difference between consistent and inconsistent trials. While both performance in consistent and inconsistent trials improved with higher scores, the performance in consistent trials grew faster (increasing from ",
              round(inv_logit_scaled(x4a$hypothesis$Estimate+ x4aa$hypothesis$Estimate),2)*100, "% accuracy with a chance level TASIT score (40) to ",
              round(inv_logit_scaled(x4a$hypothesis$Estimate + x4b$hypothesis$Estimate),2)*100, "% with a perfect score vs. an increase from ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x3aa$hypothesis$Estimate),2)*100, "% to ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x3b$hypothesis$Estimate),2)*100, "%).")

cat(text)
```

## Adjusting for IQ
```{r}
### TASIT and IQ

IQTAS_Accuracy_mo_f <- bf(
  Accuracy ~ 0 + Self : Consistency + 
    mo(EsToIQ) : Self : Consistency + 
    mo(TotAB) : Self : Consistency + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)


prior_IQTAS_Accuracy_mo <- c(
  prior(normal(1, 1), class=b),
  prior(normal(0, .3), class=b, coef=SelfOther:Consistency0:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moEsToIQ),
  prior(normal(0, .3), class=b, coef=SelfOther:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moTotAB),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

IQTAS_Accuracy_mo_m1 <- brmsFun(
  Name = "IQTAS_Accuracy_mo_m1",
  formula = IQTAS_Accuracy_mo_f,
  data = subset(dT, Group=="Schizophrenia"),
  family = bernoulli,
  prior = prior_IQTAS_Accuracy_mo
)

IQTAS_Accuracy_mo_m1_ctrl <- brmsFun(
  Name = "IQTAS_Accuracy_mo_m1_ctrl",
  formula = IQTAS_Accuracy_mo_f,
  data = subset(dT, Group=="Control"),
  family = bernoulli,
  prior = prior_IQTAS_Accuracy_mo
)
hypothesis(IQTAS_Accuracy_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  > 0", class = "bsp")
hypothesis(IQTAS_Accuracy_mo_m1_ctrl, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  > 0", class = "bsp")


IQTAS_RT_mo_f <- bf(
  RT ~ 0 + Self : Consistency + mo(EsToIQ) : Self : Consistency + 
    mo(TotAB) : Self : Consistency + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)
prior_IQTAS_RT_mo <- c(
  prior(normal(6,.5), class=b),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moEsToIQ),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moTotAB),
  prior(lkj(3), class=cor),
  prior(normal(0,.1), class=sd),
  prior(normal(0,.3), class=sigma)
)

IQTAS_RT_mo_m1 <- brmsFun(
  Name = "IQTAS_RT_mo_m1",
  formula = IQTAS_RT_mo_f,
  data = subset(dT,Accuracy == "1" & Group=="Schizophrenia"),
  family = shifted_lognormal,
  prior = prior_IQTAS_RT_mo
)

x1 <- hypothesis(IQTAS_Accuracy_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  > 0", class = "bsp")
x2 <- hypothesis(IQTAS_RT_mo_m1, "80 * SelfSelf:Consistency1:moTotAB  - 80 *  SelfSelf:Consistency0:moTotAB  < 0", class="bsp")

x3a <- hypothesis(IQTAS_Accuracy_mo_m1, "SelfSelf:Consistency0 > 0")
x3aa <- hypothesis(IQTAS_Accuracy_mo_m1, "40 * SelfSelf:Consistency0:moTotAB > 0", class="bsp")
x3b <- hypothesis(IQTAS_Accuracy_mo_m1, "80 * SelfSelf:Consistency0:moTotAB > 0", class="bsp")
x4a <- hypothesis(IQTAS_Accuracy_mo_m1, "SelfSelf:Consistency1 > 0")
x4aa <- hypothesis(IQTAS_Accuracy_mo_m1, "40 * SelfSelf:Consistency1:moTotAB > 0", class="bsp")
x4b <- hypothesis(IQTAS_Accuracy_mo_m1, "80 * SelfSelf:Consistency1:moTotAB > 0", class="bsp")

hypothesis(IQTAS_Accuracy_mo_m1, "90 * SelfSelf:Consistency0:moEsToIQ > 0", class="bsp")
hypothesis(IQTAS_Accuracy_mo_m1, "90 * SelfSelf:Consistency1:moEsToIQ > 0", class="bsp")

text <- paste0("TASIT total score was credibly related to the altercentric intrusion effect for accuracy (",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), "), but not RT performance (",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ".",
              "In particular, the higher the TASIT score, the larger the difference between consistent and inconsistent trials. While both performance in consistent and inconsistent trials improved with higher scores, the performance in consistent trials grew faster (increasing from ",
              round(inv_logit_scaled(x4a$hypothesis$Estimate+ x4aa$hypothesis$Estimate + 5.71),2)*100, "% accuracy with a chance level TASIT score (40) to ",
              round(inv_logit_scaled(x4a$hypothesis$Estimate + x4b$hypothesis$Estimate + 5.71),2)*100, "% with a perfect score vs. an increase from ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x3aa$hypothesis$Estimate + 4.53),2)*100, "% to ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x3b$hypothesis$Estimate + 4.53),2)*100, "%).")

cat(text)
```


## Animated Triangles (ATT)

```{r ATT}

ATT_Accuracy_mo_f <- bf(
  Accuracy ~ 0 + Self : Consistency + Self : Consistency : mo(ATT) + (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorATT_Accuracy_mo <- c(
  prior(normal(1,1), class=b),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moATT),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moATT),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moATT),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moATT),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

ATT_Accuracy_mo_m <- brmsFun(
  Name = "ATT_Accuracy_mo_m",
  formula = ATT_Accuracy_mo_f,
  data = subset(dT, Group=="Schizophrenia"),
  family = bernoulli,
  prior = priorATT_Accuracy_mo
)

# Is there an interaction?
hypothesis(ATT_Accuracy_mo_m, "16 * SelfSelf:Consistency0:moATT - 16 *  SelfSelf:Consistency1:moATT < 0", class="bsp")

# Calculating allocentric effect with random intentional score of 0
hypothesis(ATT_Accuracy_mo_m, "SelfSelf:Consistency0 > 0")
hypothesis(ATT_Accuracy_mo_m, "SelfSelf:Consistency0:moATT > 0", class = "bsp")
hypothesis(ATT_Accuracy_mo_m, "SelfSelf:Consistency1 > 0")
hypothesis(ATT_Accuracy_mo_m, "SelfSelf:Consistency1:moATT > 0", class = "bsp")
inv_logit_scaled(0.38 + 11 * 0.21) - inv_logit_scaled(-0.07 + 11 * 0.23)
inv_logit_scaled(0.38 + 24 * 0.21) - inv_logit_scaled(-0.07 + 24 * 0.23)
### 

ATT_RT_mo_f <- bf(
  RT ~ 0 + Self : Consistency + Self : Consistency : mo(ATT) + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorATT_RT_mo <- c(
  prior(normal(6,.5), class=b),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moATT),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moATT),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moATT),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moATT),
  prior(lkj(3), class=cor),
  prior(normal(0,.1), class=sd),
  prior(normal(0,.3), class=sigma)
)

ATT_RT_mo_m <- brmsFun(
  Name = "ATT_RT_mo_m",
  formula = ATT_RT_mo_f,
  data = subset(dT,Accuracy == "1" & Group=="Schizophrenia"),
  family = shifted_lognormal,
  prior = priorATT_RT_mo
)

Preds <- predict(ATT_Accuracy_mo_m)
dT$Preds <- NULL
dT$Preds[dT$Group=="Schizophrenia"] <- Preds[,1]

ggplot(dT, aes(ATT,Preds, color=Consistency)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  facet_wrap(.~Self) +
  theme_classic()

PredsRT <- predict(ATT_RT_mo_m)
dT$PredsRT <- NULL
dT$PredsRT[dT$Group=="Schizophrenia" & dT$Accuracy==1] <- PredsRT[,1]

ggplot(dT, aes(ATT,PredsRT, color=Consistency)) + 
  geom_point() + 
  geom_smooth(method=lm, formula=y~poly(x,2)) +
  facet_wrap(.~Self) +
  theme_classic()

# Is there an interaction?
hypothesis(ATT_RT_mo_m, "16 * SelfSelf:Consistency0:moATT - 16 * SelfSelf:Consistency1:moATT < 0", class="bsp")

# Calculating allocentric effect with random intentional score of 0
hypothesis(ATT_RT_mo_m, "SelfSelf:Consistency0 > 0")
hypothesis(ATT_RT_mo_m, "SelfSelf:Consistency1 > 0")
exp(6.54) - exp(6.72)


hypothesis(ATT_RT_mo_m, "16 * SelfSelf:Consistency0:moATT > 0", class="bsp")
hypothesis(ATT_RT_mo_m, "16 * SelfSelf:Consistency1:moATT > 0", class="bsp")
exp(6.54-0.24) - exp(6.78-0.13)

# Producing text for the manuscript

x1 <- hypothesis(ATT_Accuracy_mo_m, "16 * SelfSelf:Consistency0:moATT - 16 *  SelfSelf:Consistency1:moATT < 0", class="bsp")

x2 <- hypothesis(ATT_RT_mo_m, "16 * SelfSelf:Consistency0:moATT - 16 *  SelfSelf:Consistency1:moATT < 0", class="bsp")

text <- paste0("Accuracy in the animated triangles task (ATT) was not credibly related to the altercentric intrusion effect for accuracy (",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), ") nor was RT (",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ".")

cat(text)

```

## Relevant Psychotic Symptoms

```{r Symptoms}
SymptomsAccuracy_mo_f <- bf(
  Accuracy ~ 0 + Self : Consistency+ Self : Consistency : mo(Relevant) + 
    Self : Consistency : mo(Irrelevant) + 
    Self : Consistency : mo(SANS) + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorSymptomsAccuracy_mo <- c(
  prior(normal(1,1), class=b),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moRelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moRelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moRelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moRelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moSANS),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moSANS),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moSANS),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moSANS),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

SymptomsAccuracy_mo_m <- brmsFun(
  Name = "SymptomsAccuracy_mo_m",
  formula = SymptomsAccuracy_mo_f,
  data = subset(dT, Group=="Schizophrenia"),
  family = bernoulli,
  prior = priorSymptomsAccuracy_mo
)


Preds <- predict(SymptomsAccuracy_mo_m)
dT$Preds <- NULL
dT$Preds[dT$Group=="Schizophrenia"] <- Preds[,1]

ggplot(dT, aes(Relevant,Preds, color=Consistency)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  facet_wrap(.~Self) +
  theme_classic()
mean(dT$Preds[dT$Group=="Schizophrenia" & dT$Self=="Self" & dT$Relevant==3], na.rm=T)

# Is there an interaction?
hypothesis(SymptomsAccuracy_mo_m, 
           c("50 * SelfSelf:Consistency0:moRelevant - 50 *  SelfSelf:Consistency1:moRelevant < 0", "35 * SelfSelf:Consistency0:moIrrelevant - 35 *  SelfSelf:Consistency1:moIrrelevant < 0",  "20 * SelfSelf:Consistency0:moSANS - 20 *  SelfSelf:Consistency1:moSANS < 0"), class="bsp")

# Calculating allocentric effect with random intentional score of 0
hypothesis(SymptomsAccuracy_mo_m, c("SelfSelf:Consistency0 > 0","SelfSelf:Consistency1 > 0"))
inv_logit_scaled(2.28) - inv_logit_scaled(2.05)

hypothesis(SymptomsAccuracy_mo_m, 
           c("50 * SelfSelf:Consistency0:moRelevant > 0",  "50 * SelfSelf:Consistency1:moRelevant > 0"), class="bsp")
inv_logit_scaled(2.26+2.29) - inv_logit_scaled(2.04-3.73)

### 
SymptomsRT_mo_f <- bf(
  RT ~ 0 +  Self : Consistency + 
    Self : Consistency : mo(Relevant) + 
    Self : Consistency : mo(Irrelevant) + 
    Self : Consistency : mo(SANS) + 
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorSymptomsRTmo <- c(
  prior(normal(6,.5), class=b),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moIrrelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moRelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moRelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moRelevant),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moRelevant),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency0:moSANS),
  prior(normal(0,1), class=b, coef=SelfOther:Consistency1:moSANS),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency0:moSANS),
  prior(normal(0,1), class=b, coef=SelfSelf:Consistency1:moSANS),
  prior(lkj(3), class=cor),
  prior(normal(0,.1), class=sd),
  prior(normal(0,.3), class=sigma)
)

SymptomsRT_mo_m <- brmsFun(
  Name = "SymptomsRT_mo_m",
  formula = SymptomsRT_mo_f,
  data = subset(dT,Accuracy == "1" & Group=="Schizophrenia"),
  family = shifted_lognormal,
  prior = priorSymptomsRTmo
)

# Is there an interaction?
hypothesis(SymptomsRT_mo_m, "50 * SelfSelf:Consistency0:moRelevant - 50 *  SelfSelf:Consistency1:moRelevant > 0", class="bsp")
hypothesis(SymptomsRT_mo_m, "35 * SelfSelf:Consistency0:moIrrelevant - 35 *  SelfSelf:Consistency1:moIrrelevant > 0", class="bsp")
hypothesis(SymptomsRT_mo_m, "20 * SelfSelf:Consistency0:moSANS - 20 *  SelfSelf:Consistency1:moSANS > 0", class="bsp")

#  = -0.08, SE = 0.31, CI (-0.60, 0.42), ER = 0.65, credibility = 0.40 = -0.05, SE = 0.25, CI (-0.45, 0.37), ER = 0.71, credibility = 0.42).

###
# 

x1 <- hypothesis(SymptomsAccuracy_mo_m, "50 * SelfSelf:Consistency0:moRelevant - 50 *  SelfSelf:Consistency1:moRelevant < 0", class="bsp")
x2 <- hypothesis(SymptomsRT_mo_m, "50 * SelfSelf:Consistency0:moRelevant - 50 *  SelfSelf:Consistency1:moRelevant < 0", class="bsp")

x3a <- hypothesis(SymptomsAccuracy_mo_m, "SelfSelf:Consistency1 > 0")
x3b <- hypothesis(SymptomsAccuracy_mo_m, "SelfSelf:Consistency0 > 0")
x4a <- hypothesis(SymptomsAccuracy_mo_m, "50 * SelfSelf:Consistency1:moRelevant > 0", class="bsp")
x4b <- hypothesis(SymptomsAccuracy_mo_m, "50 * SelfSelf:Consistency0:moRelevant > 0", class="bsp")

x5 <- hypothesis(SymptomsAccuracy_mo_m, "35 * SelfSelf:Consistency1:moIrrelevant > 0", class="bsp")
x5 <- hypothesis(SymptomsAccuracy_mo_m, "35 * SelfSelf:Consistency0:moIrrelevant - 35 *  SelfSelf:Consistency1:moIrrelevant < 0", class="bsp")
x6<- hypothesis(SymptomsAccuracy_mo_m, "20 * SelfSelf:Consistency0:moSANS - 20 *  SelfSelf:Consistency1:moSANS < 0", class="bsp")

x7 <- hypothesis(SymptomsRT_mo_m, "35 * SelfSelf:Consistency0:moIrrelevant - 35 *  SelfSelf:Consistency1:moIrrelevant > 0", class="bsp")
x8 <- hypothesis(SymptomsRT_mo_m, "20 * SelfSelf:Consistency0:moSANS - 20 *  SelfSelf:Consistency1:moSANS > 0", class="bsp")

text <- paste0("Relevant psychotic symptoms were credibly related to the altercentric intrusion effect for accuracy (",
              "β = ", round(x1$hypothesis$Estimate, 2),
              ", SE = ", round(x1$hypothesis$Est.Error, 2),
              ", CI (", round(x1$hypothesis$CI.Lower, 2), ", ", round(x1$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x1$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x1$hypothesis$Post.Prob, 2), "but not for RT (",
              "β = ", round(x2$hypothesis$Estimate, 2),
              ", SE = ", round(x2$hypothesis$Est.Error, 2),
              ", CI (", round(x2$hypothesis$CI.Lower, 2), ", ", round(x2$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x2$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x2$hypothesis$Post.Prob, 2), ".",
              "). Altercentric intrusion grew with increased psychotic symptoms (from a difference in accuracy of",
              round(inv_logit_scaled(x3a$hypothesis$Estimate) - inv_logit_scaled(x3b$hypothesis$Estimate),2)*100,
              "% with 0 score, to a difference of ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x4a$hypothesis$Estimate) - inv_logit_scaled(x3b$hypothesis$Estimate + x4b$hypothesis$Estimate),2)*100,
              "% with full score). ",
              "Interestingly, performance on consistent trials increased from ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate),2)*100, "% with 0 score to ",
              round(inv_logit_scaled(x3a$hypothesis$Estimate + x4a$hypothesis$Estimate), 2)*100, "% with full score, ",
              "while performance on inconsistent trials decreased from ",
              round(inv_logit_scaled(x3b$hypothesis$Estimate),2)*100, "% with 0 score to ",
              round(inv_logit_scaled(x3b$hypothesis$Estimate + x4b$hypothesis$Estimate),2)*100, "% with full score. ",
              "Control predictors, i.e. severity of unrelated psychotic symptoms or negative symptoms were not credibly related to altercentric intrusion for Accuracy (unrelated psychotic symptoms: ",
              "β = ", round(x5$hypothesis$Estimate, 2),
              ", SE = ", round(x5$hypothesis$Est.Error, 2),
              ", CI (", round(x5$hypothesis$CI.Lower, 2), ", ", round(x5$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x5$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x5$hypothesis$Post.Prob, 2), "; SANS: ",
              "β = ", round(x6$hypothesis$Estimate, 2),
              ", SE = ", round(x6$hypothesis$Est.Error, 2),
              ", CI (", round(x6$hypothesis$CI.Lower, 2), ", ", round(x6$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x6$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x6$hypothesis$Post.Prob, 2),
              ") nor for RT (Unrelated psychotic symptoms: ",
              "β = ", round(x7$hypothesis$Estimate, 2),
              ", SE = ", round(x7$hypothesis$Est.Error, 2),
              ", CI (", round(x7$hypothesis$CI.Lower, 2), ", ", round(x7$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x7$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x7$hypothesis$Post.Prob, 2),
              "; SANS: ",
              "β = ", round(x8$hypothesis$Estimate, 2),
              ", SE = ", round(x8$hypothesis$Est.Error, 2),
              ", CI (", round(x8$hypothesis$CI.Lower, 2), ", ", round(x8$hypothesis$CI.Upper, 2), ")",
              ", ER = ", round(x8$hypothesis$Evid.Ratio, 2),
              ", credibility = ", round(x8$hypothesis$Post.Prob, 2), "."
)

cat(text)

```

## Control analysis: Symptoms and TASIT
```{r}
TasitSymptomsAccuracy_mo_f <- bf(
  Accuracy ~ 0 + Self : Consistency + 
    mo(TotAB) : Self : Consistency + 
    mo(Relevant) : Self : Consistency +
    mo(Irrelevant) : Self : Consistency +
    mo(SANS) : Self : Consistency +
    (0 + Self : Consistency | MatchID) + (1 | Picture)
)

priorTasitAccuracy_mo <- c(
  prior(normal(1, 1), class=b),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moTotAB),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moRelevant),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moRelevant),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moRelevant),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moRelevant),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moIrrelevant),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moIrrelevant),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moIrrelevant),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moIrrelevant),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency0:moSANS),
  prior(normal(0,.3), class=b, coef=SelfOther:Consistency1:moSANS),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency0:moSANS),
  prior(normal(0,.3), class=b, coef=SelfSelf:Consistency1:moSANS),
  prior(lkj(5), class=cor),
  prior(normal(0,.1), class=sd)
)

TasitSymptomsAccuracy_mo_m1 <- brmsFun(
  Name = "TasitSymptomsAccuracy_mo_m1",
  formula = TasitSymptomsAccuracy_mo_f,
  data = subset(dT, Group=="Schizophrenia"),
  family = bernoulli,
  prior = priorTasitAccuracy_mo
)

hypothesis(TasitSymptomsAccuracy_mo_m1, c(
  "80* SelfSelf:Consistency1:moTotAB > 80* SelfSelf:Consistency0:moTotAB",
  "50 * SelfSelf:Consistency1:moRelevant > 50 * SelfSelf:Consistency0:moRelevant",
  "SelfSelf:Consistency1:moIrrelevant < SelfSelf:Consistency0:moIrrelevant",
  "SelfSelf:Consistency1:moSANS > SelfSelf:Consistency0:moSANS"
), class="bsp")

```

## Figure 4: Altercentric intrusion and clinical / cognitive features

```{r}

# Create counterfactual data for tasit, att, and symptoms

dT_TASIT <- expand.grid(
  Self = "Self",
  Consistency = c("0","1"), 
  TotAB = sort(unique(dT$TotAB[dT$Group=="Schizophrenia"])), 
  Picture = unique(dT$Picture)) %>% 
  mutate(MatchID = 1000 + TotAB)

dT_TASIT1 <- expand.grid(
  Self = c("Self","Other"),
  Consistency = c("0","1"), 
  TotAB = sort(unique(dT$TotAB[dT$Group=="Schizophrenia"])), 
  Picture = unique(dT$Picture)) %>% 
  mutate(MatchID = 1000 + TotAB)

dT$Preds[dT$Group=="Schizophrenia"] <-fitted(
  TasitAccuracy_mo_m1)[,1]
dT$Preds1[dT$Group=="Schizophrenia"] <-predict(
  TasitAccuracy_mo_m1)[,1]
mean(dT$Preds[dT$Consistency==0 & dT$TotAB==37])
mean(dT$Preds[dT$Consistency==0 & dT$TotAB==78],na.rm=T)
mean(dT$Preds1[dT$Consistency==0 & dT$TotAB==37])
mean(dT$Preds1[dT$Consistency==0 & dT$TotAB==78],na.rm=T)

hypothesis(TasitAccuracy_mo_m1, c(
  "80*(SelfSelf:Consistency1:moTotAB - SelfSelf:Consistency0:moTotAB) > 80*(SelfOther:Consistency1:moTotAB - SelfOther:Consistency0:moTotAB)",
  "(SelfOther:Consistency1:moTotAB - SelfOther:Consistency0:moTotAB) >0 "
), class="bsp")

dT_TASIT$Preds <-fitted(
  TasitAccuracy_mo_m1, 
  newdata = dT_TASIT, 
  allow_new_levels = T)[,1]

dT_TASIT1$Preds <-fitted(
  TasitAccuracy_mo_m1, 
  newdata = dT_TASIT1, 
  allow_new_levels = T)[,1]

mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Self" & dT_TASIT1$Consistency==1 & dT_TASIT1$TotAB==37]) - mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Self" & dT_TASIT1$Consistency==0 & dT_TASIT1$TotAB==37])

mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Self" & dT_TASIT1$Consistency==1 & dT_TASIT1$TotAB==78]) - mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Self" & dT_TASIT1$Consistency==0 & dT_TASIT1$TotAB==78])


mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Other" & dT_TASIT1$Consistency==1 & dT_TASIT1$TotAB==37]) - mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Other" & dT_TASIT1$Consistency==0 & dT_TASIT1$TotAB==37])
mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Other" & dT_TASIT1$Consistency==1 & dT_TASIT1$TotAB==78]) - mean(dT_TASIT1$Preds[dT_TASIT1$Self=="Other" & dT_TASIT1$Consistency==0 & dT_TASIT1$TotAB==78])



dT_ATT <- expand.grid(
  Self = "Self",
  Consistency = c("0","1"), 
  ATT = sort(unique(as.numeric(as.character(dT$ATT[dT$Group=="Schizophrenia"])))),
  Picture = unique(dT$Picture)) %>% 
  mutate(MatchID = 1000 + ATT)

dT_ATT$Preds <-fitted(
  ATT_Accuracy_mo_m, 
  newdata = dT_ATT, 
  allow_new_levels = T)[,1]

dT_Psyc <- expand.grid(
  Self = "Self",
  Consistency = c("0","1"), 
  SANS = 1, 
  Irrelevant = 1, 
  Relevant = c(4, 8, 10, 12, 13, 15, 16, 26, 27, 29, 30, 31, 35, 36, 38, 45, 47), 
  Picture = unique(dT$Picture)) %>% 
  mutate(MatchID = 1000 + Relevant)

dT_Psyc$Preds <-fitted(
  SymptomsAccuracy_mo_m, 
  newdata = dT_Psyc, 
  allow_new_levels = T)[,1]

mean(dT_TASIT$Preds[dT_TASIT$Consistency==0 & dT_TASIT$TotAB==37])
mean(dT_TASIT$Preds[dT_TASIT$Consistency==0 & dT_TASIT$TotAB==78])
mean(dT_TASIT$Preds[dT_TASIT$Consistency==1 & dT_TASIT$TotAB==37])
mean(dT_TASIT$Preds[dT_TASIT$Consistency==1 & dT_TASIT$TotAB==78])


mean(dT_Psyc$Preds[dT_Psyc$Consistency==0 & dT_Psyc$Relevant==4])
mean(dT_Psyc$Preds[dT_Psyc$Consistency==0 & dT_Psyc$Relevant==47])
mean(dT_Psyc$Preds[dT_Psyc$Consistency==1 & dT_Psyc$Relevant==4])
mean(dT_Psyc$Preds[dT_Psyc$Consistency==1 & dT_Psyc$Relevant==47])

pTas <- ggplot(dT_TASIT, aes(TotAB, Preds, color = Consistency)) + 
  geom_point() + 
  #geom_smooth() +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme_classic() +
  theme(legend.position = "none") +
  ylim(0.3, 1) +
  xlab("TASIT score") +
  ylab("Predicted accuracy")

pATT <- ggplot(dT_ATT, aes(ATT, Preds, color = Consistency)) + 
  geom_point() + 
  #geom_smooth() +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme_classic() +
  theme(legend.position = "none") +
  ylim(0.3, 1) +
  xlab("ATT score") +
  ylab("Predicted accuracy")

pPsyc <- ggplot(dT_Psyc, aes(Relevant, Preds, color = Consistency)) + 
  geom_point() + 
  #geom_smooth() +
  scale_color_brewer(palette="Dark2") +
  scale_fill_brewer(palette="Dark2") +
  theme_classic() +
  theme(legend.position = "none") +
  ylim(0.3, 1) +
  xlab("Relevant psychotic symptoms score") +
  ylab("Predicted accuracy")

pTas + pATT + pPsyc
```
